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Understanding the patterns and causes of phenotypic divergence is a central goal in evolutionary biology. Much work has 
shown that mRNA abundance is highly variable between closely related species. However, the extent and mechanisms of 
post-transcriptional gene regulatory evolution are largely unknown. Here we used ribosome profiling to compare 
transcript abundance and translation efficiency in two closely related yeast species [S. cerevisiae and S. paradoxus). By 
comparing translation regulatory divergence to interspecies differences in mRNA sequence features, we show that dif- 
ferences in transcript leaders and codon bias substantially contribute to divergent translation. Globally, we find that 
translation regulatory divergence often buffers species differences in mRNA abundance, such that ribosome occupancy is 
more conserved than transcript abundance. We used allele-specific ribosome profiling in interspecies hybrids to compare 
the relative contributions of cis- and fra/75-reguIatory divergence to species differences in mRNA abundance and translation 
efficiency. The mode of gene regulatory divergence differs for these processes, as fra/w-regulatory changes play a greater 
role in divergent mRNA abundance than in divergent translation efficiency. Strikingly, most genes with aberrant tran- 
script abundance in Fl hybrids (either over- or underexpressed compared to both parent species) did not exhibit aberrant 
ribosome occupancy. Our results show that interspecies differences in translation contribute substantially to the evolution 
of gene expression. Compensatory differences in transcript abundance and translation efficiency may increase the 
robustness of gene regulation. 



[Supplemental material is available for this article.] 

The relationship between genotype and phenotype is a long- 
standing biological puzzle. Genetic mutations can affect pheno- 
types both by changing the coding sequence of proteins and by 
changing the regulation of gene expression. As much as 40 years 
ago, it became clear that differences in protein-coding sequences 
are too infrequent to explain differences in phenotype (King 
and Wilson 1975). Consequently evolution of gene regulation has 
been suggested as a major source of phenotypic diversity (Carroll 
2005; Rockman and Kruglyak 2006). 

With the advent of DNA microarrays in the 1990s, many re- 
searchers began examining gene expression differences between 
species. mRNA abundance levels were found to vary greatly both 
among and between biological species, including yeast (Brem et al. 
2002; Bullard et al. 2010; Tirosh et al. 2011), fruit flies (Ranz et al. 
2003; Rifkin et al. 2003), mice (Schadt et al. 2003), and humans 
(Pickrell et al. 2010). In general, these studies report that poly- 
morphic mRNA abundance is present in at least 25% of genes in 
a species, and mRNA levels are even more divergent between 
closely related species. However, a few recent studies suggest that 
protein abundance is more conserved than mRNA abundance 
across all domains of life, including E. coll, S. cerevisiae, and humans 
(Schrimpf et al. 2009; Laurent et al. 2010). Because these studies 
depended upon orthologous genes identifiable across vast evolu- 
tionary distance, their results could be biased toward highly 
abundant housekeeping genes. Nonetheless, these works suggest 
that post-transcriptional processes, including mRNA translation 
and protein turnover, could act to reduce the effects of variation in 
mRNA abundance and offset divergent gene expression. 
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Gene expression requires numerous steps, including tran- 
scription, mRNA translation, and protein turnover. Each of these 
processes is regulated through complex networks of ris-regulatory 
sequence elements and trans-acting factors. C/s-regulatory ele- 
ments (e.g., transcription factor binding sites and miRNA target 
sites) are localized to the alleles they regulate and thus affect gene 
expression in an allele-specific manner. In contrast, trans-acting 
factors (e.g., transcription factors, RNA binding proteins, and 
miRNAs) can affect the expression of both alleles in a diploid cell. 
Consequently, the roles of these network components can be dis- 
sected by comparing allele-specific gene expression in Fl hybrids. 
In Fl hybrids, two parental alleles are subjected to identical trans- 
regulatory environments. Thus, differences in gene expression 
from the two alleles reflect divergent activity of their c/s-regulatory 
sequences. Trans-regulatory divergence can then be inferred by 
comparing interspecific expression differences with allele-specific 
expression differences in Fl hybrids (Wittkopp et al. 2004; Tirosh 
et al. 2009). This hybrid approach has revealed that most intra- 
species variation in mRNA abundance results from trans-regulatory 
differences, but c/s-regulatory differences accumulate over evolu- 
tionary time (Lemos et al. 2008; Wittkopp et al. 2008a; Emerson 
et al. 2010). 

Divergence of gene regulatory networks has also been impli- 
cated in inheritance of gene expression levels. Like any phenotype, 
mRNA abundance can be inherited either additively or non- 
additively. One extreme form of nonadditive inheritance includes 
cases in which genes are misexpressed in Fl hybrids compared 
to parent species. This misexpression includes both over- and 
underdominant inheritance, in which total mRNA abundance in 
Fl hybrids is either much higher or lower than both parent species, 
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respectively. Misexpression of mRNA abundance can affect many 
genes (Gibson et al. 2004; Ranz et al. 2004), and is associated with 
antagonistic cis- and trans-regulatory changes (Landry 2005; 
McManus et al. 2010; Schaefke et al. 2013). 

Although the evolution of mRNA abundance has been stud- 
ied extensively, divergent regulation of mRNA translation has re- 
ceived much less attention. We used ribosome profiling to com- 
pare mRNA abundance and translation efficiency in two species 
of yeast, 5. cerevisiae and 5. paradoxus. Interspecies differences in 
mRNA abundance and translation efficiency affected similar pro- 
portions of the transcriptome. Importantly, differences in trans- 
lation efficiency tend to buffer differences in mRNA abundance. 
We further compared the relative contributions of cis- and trans- 
regulatory differences in mRNA abundance and translation effi- 
ciency using Fl hybrids. Differences in both mRNA levels and 
translation efficiency are largely attributed to trans-regulatory di- 
vergence. Finally, we investigated inheritance patterns of mRNA 
abundance and ribosome occupancy in Fl hybrids. Strikingly, ri- 
bosome occupancy was much less likely to exhibit misexpression, 
especially for genes involved in cell cycle regulation. These results 
reveal the importance of translation regulation in evolution, pro- 
vide insights into the genetic and molecular mechanisms of 
translation regulatory divergence, and indicate that translation 
regulatory networks often act as a buffer to gene regulatory di- 
vergence in yeast. 

Results 

Measuring translation efficiency with ribosome profiling 

To investigate translation regulatory divergence, we used ribosome 
profiling in diploid S. cerevisiae, S. paradoxus, and their Fl hybrid 
(Fig. 1A). This technique provides estimates of the relative number 
of ribosomes translating mRNA from each gene (ribosome occu- 
pancy). Compared with RNA-seq or microarray measurements of 
mRNA abundance, ribosome occupancy is a better proxy of protein 



abundance in yeast (Ingolia et al. 2009). This is consistent with the 
common model that yeast protein synthesis is largely regulated at 
the translation initiation step (Lackner et al. 2007; Sonenberg and 
Hinnebusch 2009; Shah et al. 2013). Ribosome-bound mRNA 
was extracted from cycloheximide-treated cells and digested with 
RNase I to produce short (28-30 nucleotides) ribosome protected 
mRNA fragments (RPFs). RPFs were then purified and used to 
generate strand-specific libraries for Illumina high-throughput 
sequencing (Ingolia et al. 2011). Strand-specific polyA-positive 
mRNA-seq libraries were made in parallel. Together, this approach 
measures both ribosome occupancy (RPF) and mRNA abundance 
(mRNA) that can be compared to estimate relative translation ef- 
ficiency, defined here as the ratio of ribosome footprints to mRNA 
fragments (RPF/mRNA) (Ingolia et al. 2009). 

We developed an analysis pipeline to compare translation 
efficiency between species and between alleles in Fl hybrids that 
leverages the high-quality genome sequences and gene annota- 
tions from each species (Scannell et al. 2011). To improve the an- 
notations of introns and translation start sites, we sequenced 
the transcriptome of S. paradoxus (Scannell et al. 2011). Our rean- 
notation of the S. paradoxus genome resulted in a working set of 
5474 orthologous protein-coding genes. Sequence reads from each 
species were aligned to their respective genomes. Reads from Fl 
hybrid samples were aligned to each species genome separately, 
and alignments were compared to identify reads that map prefer- 
entially to one allele (Fig. IB; McManus et al. 2010). This analysis 
resulted in an average of 7.4 million RPF and 6.9 million mRNA 
reads covering ORFs in each replicate (Supplemental Table SI). 

In some cases, allele-specific alignment can be biased toward 
higher quality genomes (Degner et al. 2009). However, these ef- 
fects appear to be largely absent when using two high-quality ge- 
nomes, as was the design of this study (McManus et al. 2010; Graze 
et al. 2012; Stevenson et al. 2013). To investigate the effectiveness 
of our allele-specific mapping pipeline, we combined sequence 
data from individual species to create a "mock hybrid" and com- 
pared divergence estimates from allele-specific and single-species 
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Figure 1 . Overview of allele-specific ribosome profiling to measure divergence in ribosome occupancy, mRNA abundance, and translational efficiency. 

(A) Ribosome profiling was performed on log-phase cultures of 5. cerevisiae, S. paradoxus, and their Fl hybrid. Ribosome protected fragments (RPFs) were 
purified and cloned into Illumina high-throughput sequencing libraries (left). Poly-adenylated mRNA sequencing libraries were prepared in parallel (right). 

(B) Sequence reads from each sample were aligned to both species genomes. Allele-specific reads were identified by comparing genomic alignments and 
mapped to corresponding regions of orthologous ORFs. (C) Comparison of separate and allele-specific coverage in ribosome profiling experiments. IGV 
browser tracks showing normalized coverage of RPF and mRNA sequence reads from 5. cerevisiae (blue) and 5. paradoxus (magenta) over the COX 6 gene 
(YHR051 W). Measurements of interspecies differences in read coverage are equivalent for separate species (upper) and allele-specific alignments ("mock 
hybrid," lower) (see also Supplemental Fig. SI). 
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alignment procedures. Estimates of gene expression differences 
from allele-specific and single-species alignment procedures were 
highly consistent for ribosome occupancy (RPFs; Pearson's R 2 > 
0.9 7), mRNA abundance (R 2 > 0.98) and translation efficiency 
(RPF/mRNA; R 2 > 0.96) (Fig. 1C; Supplemental Fig. SI), suggesting 
that the allele-specific analysis pipeline generates accurate estimates 
of regulatory divergence. We therefore used our allele-specific 
pipeline for both interspecies and Fl hybrid allelic comparisons. 

Interspecies differences in translation efficiency 
and mRNA abundance 

We compared interspecies differences in mRNA abundance 
and ribosome occupancy using a custom analysis pipeline (see 
Methods). Count-based tests were applied to identify statistically 
significant differences using 5% false discovery rate (FDR) and 
1.5 -fold minimum thresholds in both replicates. Our biological 
replicates were well correlated for both mRNA and RPF libraries 
(p > 0.97) (Supplemental Fig. S2). By requiring at least 20 total al- 
lele-specific mapped reads (S. cerevisiae reads + S. paradoxus reads > 
20) for both RPF and mRNA in each replicate, we identified 1308 
genes (26.9%) with differences in ribosome occupancy (RPFs), 
1739 genes (35.8%) with differences in mRNA abundance, and 
1345 genes (27.7%) with differences in translational efficiency 
(RPF/ mRNA) (Fig. 2A-C). Similar results were obtained using 
minimum cutoffs of 10, 50, and 100 reads (Supplemental Table S2). 
These results suggest that translation regulation may play as signifi- 
cant a role as mRNA abundance in the evolution of gene expression. 

Many features of mRNA transcripts influence translation ef- 
ficiency. Long transcript leaders (5' UTRs), frequent upstream 
AUGs (uAUGs), and uORF activity have all been associated with 
reduced translation efficiency (Ingolia et al. 2009; Brar et al. 2011; 
Rojas-Duran and Gilbert 2012; Waern and Snyder 2013). We used 
our RNA-seq data to annotate transcript leaders for 3014 genes in 
each species. Species differences in transcript leader length were 
negatively correlated with divergent translation efficiency such 
that the species with the longer transcript leader tended to have 
lower efficiency of translation (p = -0.27; P = 1.9 X 10" 13 ). Dif- 
ferences in the number of uAUGs were also negatively correlated 
with divergent translation efficiency, suggesting that orthologs 
with more uAUGs tend to be less efficiently translated (p = -0.27; 
P = 6.3 X 10~ 8 ). Partial correlation analysis of transcript leader 
length, controlling for uAUG frequency and vice versa, show that 
they contribute equally to divergent translation efficiency (p = 
-0.126; P= 0.00054; and p = -0.129; P = 0.00038). In comparison, 



the frequency of upstream stop codons was not correlated with 
divergent translation efficiency (P > 0.3). Thus, divergent trans- 
lation regulation is due in part to divergent transcript leader length 
and uAUG frequency. 

The importance of uAUGs suggests that upstream open 
reading frames (uORFs) may also contribute to divergent trans- 
lation efficiency. Once thought to be rare, recent studies have 
identified thousands of uORFs in 5. cerevisiae, many of which ap- 
pear to initiate translation at noncanonical "NTG" start codons 
(Ingolia et al. 2009; Brar et al. 2011). We searched for putative 
uORFs starting with NTG and ending with an in-frame stop codon 
in transcript leader regions from each species. RPF reads were 
mapped to uORF regions and the ratios of uORF to main ORF reads 
from each species were compared. Species differences in ribosome 
occupancy over uORFs were also negatively correlated with di- 
vergent translation efficiency (p = -0.25; P = 3.6 X 10~ 8 ), such that 
orthologs with higher uORF occupancy showed generally lower 
efficiency of main ORF translation (Supplemental Fig. S3). These 
results indicate that translation regulatory divergence can be at- 
tributed in part to interspecies differences in uORF presence and 
utilization. 

Due to redundancy in the genetic code, codon usage bias can 
also contribute to regulation of translation efficiency (Plotkin and 
Kudla 201 1). A previous comparison of nine yeast species revealed 
divergent codon usage bias for hundreds of orthologous ORFs 
(Man and Pilpel 2007). We calculated codon bias using the tRNA 
Adaptation Index (tAI) (dos Reis et al. 2004) to compare species 
differences in codon bias to translation regulatory differences. 
Codon bias was positively correlated to mRNA abundance, ribo- 
some occupancy, and translation efficiency (p > 0.52; P < 2 X 
10~ 16 ), such that highly expressed genes are biased toward using 
more abundant tRNA in both species. In contrast, interspecies 
differences in codon bias were weakly correlated with divergent 
translation efficiency (p = 0.11; P = 4.706 X 10" 15 ). 

Translation regulatory divergence buffers species differences 
in mRNA abundance 

Genes with interspecies differences in mRNA abundance were 
more than twice as likely to have divergent translation efficiency 
compared to genes with conserved mRNA abundance (Fisher's 
exact test [FET], P < 10~ 16 ). The co-occurrence of divergence in 
translation efficiency could work in the same (amplifying) or op- 
posite (buffering) direction (Fig. 3A). Strikingly, we found that 
buffering effects were —5.5 times more common (N = 560) than 
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Figure 2. Comparing regulatory divergence of ribosome occupancy, mRNA abundance, and translation efficiency. (A-Q Scatter plots compare the 
normalized average number of sequence reads for 5. cerevisiae (x-axis) and 5. paradoxus (/-axis). Genes with statistically significant differences in read 
counts (FDR < 5%, minimum 1 .5-fold difference) are plotted as open circles with black edges. Translation efficiency is defined here as the number of 
ribosome protected fragment reads (RPF) divided by the number of mRNA-seq reads covering an ORF. 
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Figure 3. Translation regulatory divergence buffers interspecies differences in mRNA abundance. (A) 
Cartoon depicting buffering (left) and amplification (right). mRNA are shown as blue lines, and ribo- 
somes are shown as black circles. Buffered genes have divergent mRNA abundance, with less divergent 
ribosome occupancy such that protein production is more conserved. In contrast, amplified genes have 
divergent mRNA abundance and even more divergent ribosome occupancy. (B) Scatterplot comparing 
divergent translation efficiency (/-axis, RPF/mRNA) with divergent mRNA abundance (x-axis). Buffered 
and amplified genes are plotted in red and blue, respectively. The negative correlation between mRNA 
abundance and translation efficiency suggests a genome-wide trend toward buffering. (C) Example of 
RPF and RNA-seq coverage over FPK1 , a buffered gene. IGV browser tracks showing normalized cov- 
erage for 5. cerevisiae (blue) and 5. paradoxus (red). 



amplifying effects (N = 101; FET, P < 10 ). Indeed there was 
a significant negative correlation between mRNA abundance 
and translation efficiency genome-wide (p = -0.413, P < 10~ 16 ) 
(Fig. 3B,C). This tendency toward buffering would reduce di- 
vergent expression overall at the protein level. 

We investigated several features of buffered and amplified 
genes. Highly expressed genes often have housekeeping functions 
and therefore might be buffered more frequently than genes with 
low mRNA abundance. However the mRNA abundance of buffered 
genes was either average (in S. cerevisiae) or slightly lower than 
average (S. paradoxus) (Supplemental Fig. S4). We also compared 
the biological functions of buffered and amplified genes. Gene 
ontology (GO) analysis (Supplemental Table S3) revealed that 
buffered genes were enriched in processes involved in cellular 
communication and catabolism (P < 0.01), whereas amplified 
genes, were enriched in functions related to more specific meta- 
bolic processes, including amino acid metabolism, sulfur assimi- 
lation (P < 0.01), iron import, and redox reactions (P < 0.004). The 
enrichment of these categories is consistent with species/strain 
differences in metabolic strategies, as S. paradoxus metabolism uses 
aerobic respiration more intensely, whereas S. cerevisiae has more 
auxotrophies to amino acid biosynthesis. 

Cis- and fra/w-contri buttons to translation regulatory 
divergence 

To gain further insight into the molecular and genetic changes 
underlying translation regulatory divergence, we performed allele- 
specific ribosome profiling on an Fl hybrid strain made by crossing 
our strains of S. cerevisiae and S. paradoxus. In Fl hybrids, differ- 
ences in allele-specific translation efficiency reflect divergence of 
their c/s-regulatory sequences. The effects of trans-regulatory di- 
vergence can then be inferred by comparing interspecies expres- 



sion differences with Fl hybrid allele- 
specific expression (Wittkopp et al. 2004; 
Tirosh et al. 2009; Bullard et al. 2010; 
McManus et al. 2010). This approach is 
feasible with our strains, as S. cerevisiae 
and S. paradoxus are sufficiently divergent 
to allow proper allele assignment for 71% 
of sequence reads (Supplemental Table SI). 
Figure 1C shows an example comparison of 
sequence coverage with separate species 
and allele-specific read assignments. 

At each level of gene regulatory 
control — transcript abundance, ribosome 
occupancy, and translational efficiency — 
more genes were affected by differences in 
the activity of trans-regulatory factors than 
by ds-regulatory differences (Fig. 4A). This 
is consistent with previous studies of regu- 
latory divergence and shows that the 
mode of regulatory evolution is similar for 
mRNA abundance and translation effi- 
ciency. However, processes revealed that 
ds-regulatory divergence contributes more 
to differences in translation efficiency 
(36.7% cis) than to differences in mRNA 
abundance (30.2% cis; Wilcoxon rank- 
sum test [WRT], P = 3.9 X 10" 14 ) (Fig. 4B). 

We investigated the overlap of genes 
with ds-regulatory divergence in mRNA 
abundance, translation efficiency, and ribosome occupancy (Fig. 
4C). C/s-regulatory divergence in translation efficiency was found 
three times more often in genes with c/s-regulatory divergence in 
mRNA levels than compared to those without (FET, P < 2.2 X 
10~ 16 ). Coupled c/s-regulatory evolution of mRNA abundance and 
translation efficiency was biased more than threefold toward 
buffering allele-specific ribosome occupancy. As a result, 345 genes 
(55%) with c/s-regulatory differences in mRNA levels do not ex- 
hibit c/s-regulatory divergence in ribosome occupancy. We refer to 
these genes as "c/s-mRNA buffered/' whereas genes with concur- 
rent c/s-regulatory differences in mRNA abundance and ribosome 
occupancy are "c/s-mRNA penetrant." The median magnitude of 
c/s-regulatory divergence for c/s-mRNA buffered genes is signifi- 
cantly lower than that of c/s-mRNA penetrant genes (1.84-fold 
versus 2.31-fold; WRT, P < 2.2 X 10" 16 ). This suggests that small 
magnitude c/s-regulatory differences in mRNA levels are less likely 
to penetrate into the phenotype of protein production. 

We next examined how cis- and trans-regulatory divergence 
in mRNA abundance and translation efficiency contribute to buffered 
and amplified interspecies expression differences. Compared to buff- 
ered genes, amplified genes were more likely to exhibit c/s-regulatory 
divergence of both mRNA abundance (FET, 1.6-fold, P = 0.018) and 
translation efficiency (FET, 2.3-fold, P = 5.3 X 10" 5 ). Amplified 
genes also had a higher percentage of c/s-regulatory divergence 
than buffered genes, both for mRNA abundance (31.1% versus 
25.5%; WRT, P = 0.07) and translation efficiency (60.7% versus 
31.7%; WRT, P = 0.0003342) (Fig. 4). Thus, c/s-regulatory divergence 
contributes more to amplified genes than to buffered genes. 



Buffering of misexpression in Fl hybrids 

Previous studies have revealed aberrant mRNA abundance in 
interspecies hybrids — either lower (underexpressed) or higher 



Genome Research 425 



www.genome.org 



McManus et al. 



A 2000 




Figure 4. Contributions of c/s- and trans-regulatory divergence in mRNA abundance, ribosome oc- 
cupancy, and translation efficiency. (A) Bar plot shows the number of genes affected by significant 
regulatory divergence in c/'s-acting sequences (C) and trans-acting factors (T). Irons-acting factors 
contribute most to mRNA abundance. (B) Box plot showing the fraction of regulatory divergence at- 
tributable to differences in c/s- regulatory elements (%c/s). Ribosome occupancy and translation effi- 
ciency both have higher %cis than mRNA abundance. (C) Venn diagram showing the overlap of genes 
with c/s- regulatory divergence in mRNA abundance, ribosome occupancy, and translation efficiency. 
(D) Box plot depicting %cis for genes with buffering (Buf) and amplifying (Amp) regulatory divergence 
of mRNA abundance (left) and translation efficiency (right). Asterisks indicate results of Wilcoxon rank- 
sum test comparisons: (*) P < 0.001 ; (**) P < 0.0005. 



(overexpressed) than both parental species. This misexpression 
may be related to interspecies network incompatibilities and, po- 
tentially, to speciation (Ranz et al. 2004; Landry 2005; Moehring 
et al. 2007). However, a recent proteomics study estimated that 
only 3% of 398 analyzed genes exhibit misexpression of protein 
abundance in Fl hybrids of 5. cerevisiae and S. bay anus (Khan et al. 
2012). We analyzed the inheritance of divergent expression in 
both mRNA abundance and ribosome occupancy to test the hy- 
pothesis that mRNA misexpression could be translationally buff- 
ered. Compared to mRNA abundance, ribosome occupancy showed 
strikingly fewer genes with nonadditive inheritance modes (Fig. 5), 
with 6.3-fold fewer underexpressed genes (FET, P < 10~ 16 ) and 1.9- 
fold fewer overexpressed genes (FET, P = 0.0003194). This result 
indicates that divergent translation regulation also buffers mis- 
expression of mRNA abundance, such that total protein synthesis 
in Fl hybrids is more consistent with that of parent species. 

Buffering of misexpression could reduce hybrid incom- 
patibilities. To investigate the effects of misexpression and buff- 
ering, we compared enrichment of GO functions in genes misex- 
pressed at the mRNA level and at the level of ribosome occupancy. 
Genes underexpressed at the mRNA level were functionally 
enriched in cell cycle regulation and cellular metabolic processes, 
whereas genes with overexpressed mRNA largely functioned in 
purine biosynthesis pathways (P < 0.01) (Supplemental Table S4). 
Consistent with this, genes with aberrantly high ribosome occu- 
pancy in the Fl hybrid strain were enriched for functions in inosine 
monophosphate biosynthesis, the precursor to purine synthesis. In 
contrast, no GO functional categories were enriched in the 34 genes 
underexpressed at the translation level. In conclusion, compensa- 
tory changes in translation efficiency appear to correct aberrant 
mRNA abundance in Fl hybrids of S. paradoxus and S. cerevisiae, 
possibly restoring the expression of cell cycle control genes. 

Discussion 

The evolution of translation efficiency 

Although undoubtedly important, evolutionary changes in mRNA 
abundance may not always reflect interspecies differences in gene 



expression. Indeed, prior work has high- 
lighted a disconnect between variation in 
mRNA and protein abundance within 
yeast (Foss et al. 2007), mice (Ghazalpour 
et al. 201 1), and humans (Wu et al. 2013). 
Together, these studies suggest that vari- 
ation in mRNA translation and pro- 
tein turnover contribute to polymorphic 
gene expression. Our results indicate that 
translation regulation plays a substantial 
role in the evolution of gene expression. 
Roughly one-quarter of all expressed genes 
exhibited divergent translation efficiency, 
showing that divergence in translation 
regulation and mRNA abundance affect 
similar numbers of genes. 

Changes in codon usage have been 
previously implicated in yeast gene reg- 
ulatory evolution (Man and Pilpel 2007). 
More recent studies have highlighted 
important translation regulatory roles of 
transcript leaders in S. cerevisiae (Ingolia 
et al. 2009; Brar et al. 2011; Rojas-Duran 
and Gilbert 2012; Waern and Snyder 2013). Our results indicate 
that interspecies differences in transcript leaders appear to play 
a larger role than codon bias. One explanation for this observation 
is that changes in codon bias are more likely to be pleiotropic, as 
changes in codon usage can alter translation elongation and dis- 
rupt cotranslational protein folding (Angov 2011; Plotkin and 
Kudla 2011). In contrast, changes in transcript leader length and 
sequence composition are likely to affect only translation initia- 
tion. Much like transcription promoter regions, transcript leaders 
provide important opportunities for evolutionary changes in gene 
expression without affecting protein-coding sequences. 

Differences in modes of regulatory evolution for mRNA 
abundance and translation efficiency 

We used allele-specific ribosome profiling to compare the contri- 
butions of cis- and trans-regulatory changes toward divergence of 
mRNA abundance and translation efficiency. We found that cis- 
regulatory differences appear to contribute more to divergence in 
translation efficiency than to interspecies differences in mRNA 
abundance. This is inconsistent with other work suggesting that 
fewer ds-acting loci contribute to polymorphic regulation of 
mRNA abundance than of protein abundance among S. cerevisiae 
strains (Foss et al. 2011). This may be due in part to a smaller 
sample size (354 genes) and less power of mass spectrometry-based 
studies for QTL mapping. In fact, recent work from the same group 
has shown that local pQTL are more common than previously 
appreciated (FW Albert, S Treusch, AH Shockley, JS Bloom, L Kruglyak, 
in prep.). Differences in divergence time in the Foss study (between 
strains) and ours (between species) could also contribute to their 
differing conclusions, as it has been shown that ds-regulatory di- 
vergence accumulates over evolutionary time (Lemos et al. 2008; 
Wittkopp et al. 2008b; Emerson et al. 2010). 



Buffering of divergent mRNA abundance 
and translation efficiency 

Our results indicate that regulatory interactions between mRNA 
abundance and translation efficiency more often reduce than 
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Figure 5. Inheritance of gene expression in F1 hybrids of 5. cerevisiae and 5. paradoxus. (A) Hypothetical patterns of gene expression in 5. paradoxus 
(red), 5. cerevisiae (blue), and F1 hybrid yeast (purple), depicting six modes of gene expression inheritance. (B) Scatterplots showing the difference 
between mRNA expression levels (left) in the F1 hybrid and 5. cerevisiae (x-axis) and 5. paradoxus (/-axis). The difference between ribosome occupancy 
levels (RPF) in the F1 hybrid and parental species is shown on the right. Bar plots show the number of genes in each inheritance category from ribosome 
occupancy data. Genes with overdominant (overexpressed) and underdominant (underexpressed) inheritance occur much less frequently when con- 
sidering ribosome occupancy, as compared to misexpression in mRNA abundance. 



amplify interspecies differences in gene expression. Buffering of 
gene expression may help explain recent work showing that pro- 
tein abundance is more conserved than mRNA abundance across 
organisms (Schrimpf et al. 2009; Laurent et al. 2010). Changes in 
translation efficiency may compensate for changes in mRNA 
abundance by balancing the numbers of translating ribosomes. 

There are two potential, nonexclusive mechanisms for buff- 
ering divergent gene regulation. First, regulatory networks con- 
trolling mRNA abundance and translation efficiency could have 
diverged through compensatory genetic mutations ("genetic 
compensation"), such that mutations in one pathway counteract 
the effects of mutations in a different pathway. Secondly, the gene 
regulatory networks may be inherently robust, such that changes 
in mRNA abundance are automatically buffered by feedback loops 
and regulatory circuitry ("network robustness"). The contributions 
of cis- and trans-regulatory changes in buffered genes may help 
differentiate between these mechanisms of buffering. For example, 
the genetic compensation model predicts that buffered genes 
would more frequently exhibit counteracting ds-regulatory changes 
in mRNA abundance and translation efficiency. In contrast, coun- 
teracting trans-acting differences are more consistent with network 
robustness. Of the 560 buffered genes in our data set, 75 (13%) 
exhibit counteracting ds-regulatory changes, whereas 254 (45%) 
have counteracting trans-acting changes. Consequently, we pro- 
pose that network robustness is likely to be more common than 
genetic compensation. 

Compared to buffered gene expression levels, amplified gene 
expression levels have fewer concordant trans-regulatory changes 
in both mRNA levels and translation efficiency. Of the 101 am- 
plified genes, 23 (23%) have amplifying ds-regulatory changes in 
both mRNA abundance and translation efficiency, whereas 24 
(24%) have amplifying trans-regulatory changes. Amplifying genes 
also had a higher percentage of ds-regulatory divergence, showing 
that ds-regulatory changes contribute more to divergent regula- 
tion of amplified genes, compared to buffered genes. Importantly, 



ds-regulatory divergence is thought to be more frequently mediated 
by positive selection than trans-regulatory divergence (Emerson 
et al. 2010; Graze et al. 2012; Schaefke et al. 2013). Thus amplifying 
divergence of mRNA abundance and translation efficiency may re- 
flect positive selection for gene expression, as has been suggested for 
combinations of ds-regulatory mutations affecting mRNA abun- 
dance alone in Drosophila (Graze et al. 2012). 

Many biological systems, including gene regulatory net- 
works, are known to exhibit robustness to genetic and environ- 
mental perturbation (Stark et al. 2005; Levy and Siegal 2008; Masel 
and Siegal 2009; MacNeil and Walhout 2011; Denby et al. 2012). 
By increasing the robustness of protein production, translational 
buffering could reduce the phenotypic impacts of variation in 
mRNA abundance. This in turn would allow increased variation in 
mRNA abundance which, if unmasked by disruption of buffering, 
could give organisms more variation upon which selection might 
act. This would create a situation identical to canalization in de- 
velopment (Waddington 1942; Rutherford and Lindquist 1998), in 
which capacitor genes allow the buildup of latent variation, which 
is only observed under extraordinary conditions. Our results sug- 
gest that translational buffering serves a similar purpose by allowing 
variation in mRNA abundance to accumulate through genetic drift 
without substantially changing protein levels. 

Earlier studies on gene regulatory evolution have documented 
widespread misexpression of mRNA abundance in interspecies hy- 
brids (Gibson et al. 2004; Ranz et al. 2004), often associated with 
counteracting cis- and trans-regulatory divergence of mRNA abun- 
dance (Landry 2005; McManus et al. 2010; Schaefke et al. 2013). 
Surprisingly, our analysis revealed translational buffering of misex- 
pressed genes. Hundreds of genes had aberrant mRNA levels in 
hybrid cells, yet few showed aberrant ribosome occupancy pat- 
terns. As with interspecies differences in mRNA abundance, over- 
and underexpressed genes could be buffered through genetic 
compensation or network robustness. It has been shown that 
protein-protein interactions are functional in Fl hybrids of 
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S. cerevisiae and S. kudriavzevii (Leducq et al. 2012), suggesting that 
many biochemical features of regulatory networks are robust in 
interspecies hybrids. Thus, network robustness is an attractive 
answer. By reducing aberrant gene expression patterns, trans- 
lational buffering of misexpression could contribute to the con- 
siderable success of yeast hybrids (Morales and Dujon 2012). 

Thus far, we have investigated translation regulatory evolu- 
tion only in the context of optimal growth conditions. Because 
many stress responses involve changes in mRNA translation, re- 
peating these experiments under different conditions could shed 
light on the role of translation in species-specific responses to 
stress. It will also be important to identify the genes required to 
maintain translational buffering, as they could act as capacitors 
for variation in mRNA abundance. Although there is much yet to 
learn, our results underscore the importance of translation in the 
evolution of gene expression, help define the molecular and ge- 
netic causes of divergence in translation efficiency, and show that 
interspecies differences in translation efficiency often act as a 
buffer to gene regulatory evolution. 

Methods 

Yeast strains 

Saccharomyces strains: S. cerevisiae (GSY83; S288C), S. paradoxus 
(GSY82; CBS432), and their Fl hybrid (GSY88) were kind gifts from 
Gavin Sherlock (Kao et al. 2010). The original GSY83 is a haploid 
strain, whereas GSY82 and GSY88 are diploid strains. For this 
study, a diploid version of GSY83 was produced by galactose- 
induced transient expression of the HO gene from plasmid 
YCP50-Gal::HO, a gift from John Woolford. Diploid GSY83 were 
identified via PCR screening. 

Ribosome profiling library preparation 

Yeast polysome extracts were prepared for ribosome profiling as 
previously described (Ingolia 2010) with slight modifications (see 
Supplemental Methods). RNA for matched RNA-seq was extracted 
from 50 |xL of cell lysate using TRIzol (Life Technologies) and 
vortexing with glass beads. The sample was then enriched for 
mRNA using magnetic oligo-DT DynaBeads (Life Technologies) 
according to the manufacturer's instructions. The resulting RNA 
was then fragmented by incubating in alkaline fragmentation 
buffer (50 mM NaCO s , pH 9.2; 1 mM EDTA) for 40 min at 95°C, 
precipitated with one volume of isopropanol and resuspended in 
20 (jlL of nuclease-free water. 

High-throughput sequencing library preparation was per- 
formed as previously described (Ingolia et al. 2011) with slight 
modifications (see Supplemental Methods). RPF and mRNA library 
cDNA templates were amplified by 12 cycles of PCR using Phusion- 
polymerase (New England Biolabs), with primers incorporating 
barcoded Illumina TruSeq library sequences. The resulting PCR 
products were then purified using PCR purification columns 
(Qiagen) according to the manufacturer's instructions. The quality 
and size of the PCR products were assessed using an Agilent 
Tapestation. Libraries were sequenced on an Illumina HiSeq 2000 
at the University of Southern California Epigenome Center. 

Sequence data alignment 

Raw sequences were trimmed of 3' adapters using custom Perl 
scripts. For RPF data, reads with lengths of 27-33 nucleotides were 
retained for further analysis as this size is most likely to represent 
footprinted fragments. For mRNA, reads with lengths of 27-40 



nucleotides were retained. Adapter trimmed reads were first 
aligned to genomes using Bowtie (Langmead et al. 2009) (version 
0.12.8) with parameters -best -chunkmbs 500. Reads that failed 
to map to genomic regions were recovered and aligned to 
splice-junction databases requiring a minimum overlap of 6 nt. 
Splice-junction read alignments were converted to SAM formatted 
genomic coordinates using a custom Perl script. 

Allele-specific read assignment was performed as previously 
described (McManus et al. 2010) with slight modifications. Briefly, 
adapter trimmed reads from separate species were concatenated 
into single "mixed parent" files. Mixed parent and hybrid se- 
quence files were aligned to both species genome and junction 
databases as described above, and alignments were compared with 
a custom Perl script to identify reads that map preferentially to one 
genome (Fig. IB). The library preparation procedure used here 
frequently adds untemplated nucleotides during reverse tran- 
scription (Brar et al. 2012). Considering this, we ignored mis- 
matches at the first two bases of the 5' end of sequence reads when 
comparing alignments between two species. Reads were assigned 
to the allele to which they map with the fewest mismatches. Esti- 
mates of interspecies differences in gene expression made through 
our allele-specific pipeline were highly conelated with those made 
by single-species alignment used in traditional RNA-seq and ribo- 
some profiling (p > 0.96 for RPF, mRNA, and translation efficiency) 
(Supplemental Fig. SI). 

Gene expression analyses 

Allele-specific reads were mapped to ORFs and uORFs and nor- 
malized by down-sampling all genes from the species or allele with 
more mapped reads (McManus et al. 2010). We also accounted for 
species differences in ORF length by down-sampling reads from 
the longer ortholog in proportion to the ratio of ortholog lengths 
(i.e., if one species' ORF was twice as long, we would divide the 
number of assigned reads by two). Tests of statistical significance 
were performed essentially as previously described (McManus et al. 
2010) using R (version 2.15.2). For analyses of expression differ- 
ences, raw P- values were corrected using the Benjamini-Hochberg 
method (Benjamini and Hochberg 1995), and significant differ- 
ences were identified using 5% FDR and 1.5-fold magnitude 
thresholds. Significant differences in mRNA abundance and ribo- 
some occupancy (RPF) were identified using binomial exact tests 
(BET), requiring both replicates to meet significance thresholds. 
Differences in translation efficiency (RPF/mRNA) were identified 
using Cochran Mantel-Haenszel tests (CMH), which directly in- 
corporate experimental replicates. CMH tests were used to evaluate 
trans-regulatory divergence by comparing the ratio of expression 
differences in parents and hybrids for mRNA abundance and ribo- 
some occupancy. The method of Altman and Bland (2003) was 
used to test for trans-regulatory divergence in translation effi- 
ciency. Supplemental Table S5 contains all read-count data, tAI 
calculations, expression ratios, and results of significance test for 
divergent expression. For full expression analysis details, see the 
Supplemental Methods. 

Genes were classified as "buffered" or "amplified" using the 
following criteria. Both buffered and amplified genes have signif- 
icant differences in mRNA abundance and translational efficiency. 
In buffered genes, the interspecies difference in ribosome occu- 
pancy (RPF) is at least 1.5 -fold lower than the interspecies differ- 
ence in mRNA abundance. In contrast, amplified genes have at 
least 1.5-fold higher differences in ribosome occupancy compared 
to differences in mRNA abundance. 

Categories of gene expression inheritance were determined 
using a custom R program as previously described (McManus et al. 
2010). Total expression in Fl hybrid cells was normalized to that of 
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parental samples and log- transformed. Genes with total expression 
more than 1.25-fold different from that of either parent species 
were considered to have nonconserved inheritance and were 
classified as having additive, dominant, underdominant, or over- 
dominant inheritance modes based on the magnitude of differ- 
ences between hybrid and parental expression. 

Additional analyses 

Gene ontology enrichment analysis was performed using the 
generic gene ontology term finder (http://go.princeton.edu), and 
GO terms were further evaluated using revigo (http://revigo.irb.hr) 
(Supek et al. 2011). The tRNA adaptation index (tAI) was calcu- 
lated using the codonR package (http://people.cryst.bbk.ac. 
uk/~fdosr01/tAI/index.html) (dos Reis et al. 2004). Genome 
Browser tracks (bedGraph) were produced for visual analysis using 
genomeCoverageBed from BEDTools. Images for all figures were pro- 
duced using the Integrative Genome Browser (Robinson et al. 2011). 

Data access 

High-throughput sequencing data have been submitted to the 
NCBI Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/ 
sra) under accession number SRP028552. A read-count table for 
mRNA and RPFs, genome-browser tracks (bedGraph), and gene 
annotations (bed) have been submitted to the NCBI Gene Ex- 
pression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) un- 
der accession number GSE52119. 
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